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A method is described for building an equation of state which <iivcs an infinite specific heat C t at 
the critical point approached from the one-phase domain. Interpolated PVT data are used on isochores 
up to 100 K at densities to the liquid triple point. In addition to an accurate representation of these 
data, the equation #ives specific heats which agree well with experimental data along near-critical 
isochores. 
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Symbols and Units 

(Subscripts c and t refer to critical and liquid triple 
points) 

A, B, C, D, density-dependent coefficients of eq (5). 

C v (p, T) , specific heat, J/mol K. 

C ( l(T) , specific heat at p = 0, J/mol K. 

C*(p, T) , interaction specific heat, J/mol K. 

d, density. 

d c , 15.59 mol/1 for eq (3), [13] . 

d u (1.0/0.026176) mol/1, [19 1. 

,/, the joule. 

/, the liter. 

mol, the gram mole, based on R in this work 

[5]. 

P, pressure, 1 atm = 0.101325 MN/m 2 . 

R, the gas constant, 0.0820597 l-atm/(mol- 

deg) [5]. 

p = d/dt, density reduced at the liquid triple point. 

cr = d/d c , density reduced at the critical point. 

T, temperature K, (NBS 1955). 

7V, 32.93 K. 

T u 13.803 K, [19]. 

T s (p) , coexistence temperature, K. 

0(p), a locus of temperatures, K. 

v, ljd, molal volume. 

co(p, T)= l-0(p)/T. 

x = T/T t , temperature reduced at the triple point. 

Z=Pv/RT, the "compressibility factor". 



1. Introduction 

The objective of this work is an accurate form of an 
equation of state to replace costly approximation 
methods in the computation of thermodynamic prop- 
erties [1, 2] 1 . As experimental specific heats C r 
increase without limit on approach to the critical point 
[3, 4], this behavior is included in the present equa- 
tion. The nonanalytic term is selected to conform with 
specific heats of parahydrogen at temperatures 
relatively far removed from the critical point. The 
general form of this isochoric equation of state can be 
related to the well-known Beattie-Bridgeman equation. 

Three particular benefits derive from the present 
selection of temperature-dependent terms in the 
equation. Coefficients of these terms (found independ- 
ently on each isochore of PVT data) are smooth func- 
tions of the density. Given the critical temperature, 
the present equation yields a critical density (and 
pressure) from PVT data in agreement with the known 
critical density of parahydrogen. Specific heats from 
the present equation agree with experimental data 
along near-critical isochores, without the use of 
complicated adjustments to the equation for the critical 
region. One may hope that present methods can be 
applied to other substances. 

Data used in this work are from [5, 6]. Excellent 
reviews on critical phenomena are available [7, 8]. 

2. The Behavior of Specific Heats 

Specific heats are related to (d 2 PldT 2 ) f > from the 
equation of state by integration along isotherms: 
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c*=-T-l\d 2 Pidr)dp/ P 2 

Jo 

C* = C v (p 9 T)-C° v (T) 



(1) 



<o( P ,T)=l-0(p)IT, 



(2) 



is the interaction specific heat, and C$(T) at p = 
refers to ideal gas states (translational and intra- 
molecular effects only). 

With any equation of state, eq (1) generally will 
yield several temperature-dependent terms for C*. 
As T^>T ( along the critical isochore, however, there 
may be a nonanalytic term which exceeds all others. 
This term may be algebraic (T— T c )~ a , [9], or it 
may be logarithmic, -log {T~T C ), [10]. Assuming 
the algebraic form, we examine some parahydrogen 
data in figure 1 at p/p c =1.05. The slope of this 
log-log plot corresponds to exponent a. For these 
data at T> 1.03 • T c the exponent a = 3/4 is not com- 
parable with a =0.05 deduced for other substances at 
T< 1.03-7; [11]. We have used the argument co 
in figure 1 because integration will yield a power of 
T not exceeding unity in the equation of state. As 
background for figure 1, our oxygen specific heats 
were represented by use of a= 1/3 together with other 
terms [12]. 

3. An Origin for Temperature on Isochores 

Consider the expression 

C?~[(T-T e )IT]- a 

at densities departing from the critical isochore D — C 
of figure 2. We must replace T c by a locus of tempera- 
tures 6(p) inside the coexistence envelope [12] 
which may serve as the origin for temperature along 
each isochore. If now we define the useful function 
of temperature along each isochore, 




-Xoq CO 



Figure 1. Interaction specific heats C* {J/mol K) for p — H 2 . 

Open circles from [6|; filled circles from [2]; barred circle from [15]. 



we have in (l/co) a pole at the critical point approached 
from the one-phase domain. 

To build an equation of state we need a formula for 
d(p). An absolute description is difficult to formulate 
under the condition ^ T s . For present work, therefore, 
I define 6 as a fraction of the saturation temperatures 
formulated in the appendix. The following is selected 
from several definitions investigated: 



0(p)=r s (p)-exp[-a-(cr-l) 2 ]. 



(3) 



As the calculated T s (p) terminates at the liquid triple 
point, we must for the present eliminate PVT data at 
higher densities. Equation (3) is arbitrary, but it 
permits us to go ahead with the formulation of an 
equation of state. It obviously deserves further in- 
vestigation. It is interesting to note that the van der 
Waals' school uses a similar locus defined by maxima 
and minima of calculated isotherms inside the 
saturation envelope [14]. 

4. Formulation of the Nonanalytic Term 

Consider the calculation of specific heats by means 
of (1) along the critical isotherm A — C — B of figure 2. 
As C v is at a maximum at point C, it is necessary that 
d 2 P/dT 2 change sign as density increases through the 
critical point. Assuming that one term alone in the 
equation of state is significant near point C, the change 
of sign may be given by a coefficient, D(p), in the 
abbreviated equation: 

P~D(p)-il,(p,T) 

where d^/dT 2 ^™asT^>0 + . 
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FIGURE 2. Density-temperature phase diagram. 

T s (p) is the liquid-vapor coexistence envelope. Point C is the critical point. 
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Consider next that PJ/T isochores at p < p c have only 
negative curvature, whereas at p > p c they have posi- 
tive curvature at low temperature, changing to negative 
curvature at high temperature [15]. To facilitate the 
construction of a wide-range equation of state, let the 
nonanalytic function have the same qualitative be- 
havior as low-density isochores, figure 3. 

Nonanalytic forms with monotonically negative 
curvature include: 



i//(p, T) = o) — wlog O), 



(4a) 



i//(p, T) - [a)-co 2 -«/(2-a)]/[l-l/(2-c0], (4b) 



^(p 9 T) = co-f 5-w 2 -log co, 8 = const. 



(4c) 



All of these forms when used in (5) below give equally 
good representation of PVT data provided we omit the 
first point of every isochore (within 1 K of saturation). 
With (4c), pressures on the near-critical isotherm at 
33 K cannot be represented acceptably. Function 
(4b) gives improving representation of PVT data as 
a->l, leading to the selection of function (4a) for 
the equation of state. 

5. The Equation of State 

Our method is to postulate temperature-dependent 
terms in the equation of state. Coefficients of these 
terms are found by least squares independently on 
each isochore of PVT data. If, in addition to an accur- 
ate representation of the data, we obtain coefficients 
which are regular (smooth) functions of density, we 
may hope to develop the extremely precise descrip- 
tions required for these coefficients. 




.8 2.0 




FIGURE 3. Two different functions i//(p, T), nonanalytic at 7 = 0. 



Equation of state (5) evolves from many exploratory 
computations with the interpolated PVT data on 
isochores [5]. Symbols are listed elsewhere. This equa- 
tion uses co(p, T) defined by (2), 6(p) defined by (3), 
and density-dependent coefficients A, B, C, D: 

(Z-l) • x/p = A + B-(S>(T) 

+ C-X(p, T)+D-V(p 9 T). (5) 

The temperature-dependent functions in (5) are 

®(T) ■ x • [1-6 -exp (-£/*)], (5b) 

X(p,T)=0(p)IT, (5c) 

*(p, T) =l-[o>-a>-loga>]. (5d) 

Function O(T) is designed for freedom in high-tem- 
perature description of the second virial coefficient 
and to give at the same time a sensible behavior as 
T—>0. The condensation of two terms into one is 
necessary to minimize the total number of terms in (5). 

Function X(p, T) was selected by trial both for 
simplicity in C(p) and for representation of PVT data. 
As we shall have D(p)=0 on the critical isochore, it 
follows that the inflection on this isochore for para- 
hydrogen is given by <£>(T) and X(p, T) acting in 
concert. Efforts to condense this behavior into a single 
term have failed. 

Function "fy(p,T) gives the nonanalytic behavior 
02^/37*2 _> oo as r->0, i.e. as co^O. Through defini- 
tion (5d) this function approaches zero at high tempera- 
ture. Its formulation was discussed above. At T>6 
we have ¥-> (O/T) 2 . 

Density-dependent coefficients of (5) at first are 
obtained by least squares independently on each iso- 
chore of PVT data. Extreme precision is required in 
the representations of these coefficients because (5) 
gives a small difference of large values. For this reason 
we proceed in steps as follows: 

(1) The original data for B(p) are seen in figure 4b. 
We select a linear form; transpose the calculated term 
B(p) - <£>(T) to the left side of (5); and redetermine 
A(p), C(p), D(p) on each isochore. Figure 4a gives 
these data for A(p). 

(2) A third order power series is selected for A(p) 
after testing in combination with other variables in (5). 
With the calculated A(p)+B(p) • &(T) on the left 
of (5) we redetermine data for C(p) and D(p). 

(3) We attack D(p) next due to our special interest 
in the root near p c . The data of figure 4d suggest a 
cubic, but a fifth order series is necessary to reproduce 
PVT data. With three calculated terms transposed to 
the left of (5) we average C(p) on each isochore, with 
results in figure 4c. 

For C(p) a seventh order series is found necessary 
by trial. All of these results are summarized as follows: 

A(p) = j^Arp l , C(p) = j?Cr P ', 
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B{p) = ^Brp\ D(p) = 2D i -p> 



6. Discussion of Results 

Interpolated data on 40 isochores are used at 
densities through 37.0 mol/1. The first point of each is 
omitted, leaving 953 PVT points for least squares 
determination of coefficients [5]. 

The hidden constants a, b, (3 in (5) may be investi- 
gated at various stages in formulating (5). Table 1 
gives a survey using the above power series for 
coefficients A, B, C, D. For each value of a, b, and f3 
we observe the mean relative deviation of 953 pressures 
in parts per ten thousand "PCT* 100" and also the 
root in D(p), "D, ROOT" in mol/1. For values a = 0.2, 
6 = 0.39, /3=10, the deviation is at a minimum 0.0203 
percent and the root at 15.54 mol/1 is in excellent 
agreement with 15.27 ±0.30 mol/1 from [16] and with 
15.59 ±0.05 from [2]. 

Table 2 gives all of the constants for equation of 
state [5]. At this critical point, 32.93 K and 15.54 mol/1, 
we obtain the following comparison with our vapor 
pressure equation [17] 









d log P\d 




P, atm 


P, MN/m 2 


log/ 1 


Equation (5) 


12.628 


1.2795 


4.960 


V.P. eq, [17] 


12.673 


1.2841 


4.833 



Table 3 gives calculated pressures (relative to the 
critical pressure) along the critical isotherm at 32.93 K. 
The small imperfection of this isotherm could arise 
for many reasons, including the need for a smaller 
critical exponent a at temperatures very near T c . 

The most exacting check on (5), however, is the com- 
parison of specific heats. We have selected data 
along near-critical isochores for this purpose. In tables 
4 and 5 the data of Younglove and Diller are diminished 
to C*. Density along the isochores is given on the first 
row. Columns from left to right give T and CJ; then C* 
from [6] and as calculated with (5); and last the differ- 
ence of these values. 

The authors ascribe an uncertainty of "several 
percent" at temperatures near saturation. The greatest 
difference, 0.46 J/mol K at 36.035 K in table 5, is 3.0 
percent of the experimental C v . It is interesting to 
make similar comparison with calculated specific heats 
from our monograph [2] in table 6 to see the increasing 
deviations as 7 1 — > T c . Tedious and costly adjustments 
were required in [2] to conform with experimental C v 
data, whereas (5) now gives acceptable values directly. 

The Beattie-Bridgeman equation [18] may be written 
in the abbreviated, isochoric form 



(Z-l) -xlp = A(p)+B(p) -x + D(p)lx 2 . 



(7) 



We note absence from (7) of a term in l/x like X(p, T) 
in (5). This might be related to the fact that C{p) in (5) 
is the only coefficient with origin near zero, figure 4c. 
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FIGURE 4. Density dependence of coefficients A, B, C, D. 
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Table 1. Survey for constants a, b, and (3 



Table 3. Calculated critical isothe 







a = 


0.10 






P 


b 


0.38 


0.39 


0.40 


0.41 


8.00 


PCT* 100 = 


3.46 


3.62 


3.94 


4.36 


8.00 


D, ROOT = 


14.41 


14.61 


14.87 


15.12 


9.00 


PCT* 100 = 


2.62 


2.88 


3.29 


3.81 


9.00 


Z>, ROOT = 


14.90 


15.14 


15.40 


15.68 


10.00 


PCT* 100 = 


2.16 


2.33 


2.69 


3.21 


10.00 


Z), ROOT = 


15.18 


15.44 


15.70 


15.98 


11.00 


PCT* 100 = 


2.57 


2.39 


2.56 


2.88 


11.00 


D. ROOT = 


15.28 


15.53 


15.79 


16.06 






a = 


0.20 









b 


0.38 


0.39 


0.40 


0.41 


8.00 


PCT* 100 = 


3.32 


3.62 


4.17 


4.85 


8.00 


Z), ROOT = 


14.52 


14.73 


14.96 


15.20 


9.00 


PCT* 100 = 


2.31 


2.63 


3.22 


3.95 


9.00 


Z), ROOT = 


15.01 


15.25 


15.50 


15.77 


10.00 


PCT* 100 = 


2.25 


2.03 


2.32 


2.88 


10.00 


D, ROOT = 


15.30 


15.54 


15.80 


16.08 


11.00 


PCT* 100 = 


3.51 


3.01 


2.67 


2.71 


11.00 


D, ROOT = 


15.41 


15.64 


15.89 


16.16 






a = 


0.30 









b 


0.38 


0.39 


0.40 


0.41 


8.00 


PCT* 100 = 


3.56 


4.09 


4.95 


6.04 


8.00 


D, ROOT = 


14.61 


14.80 


15.01 


15.24 


9.00 


PCT* 100 = 


2.50 


2.70 


3.57 


4.65 


9.00 


D, R()OT = 


15.11 


15.33 


15.57 


15.84 


10.00 


PCT* 100 = 


3.17 


2.35 


2.27 


2.88 


10.00 


D, R()OT = 


15.41 


15.64 


15.89 


16.17 


11.00 


PCT* 100 = 


5.27 


4.45 


3.72 


3.16 


11.00 


D, ROOT = 


15.54 


15.76 


15.99 


16.26 



Table 2. Constants of eq (5) for parahydrogen 



a = 0.2, 6 = 0.39 = 10 



Ao = 


-9.1083 


1777 


c 4 =- 


-4S8.5219 


0867 


A,= 


-2.3668 


5605 


c 5 = 


606.6850 


4566 


A 2 = 


5.5797 


7042 


c«=- 


-416.4729 


0930 


A,= 


2.4230 


4292 


c 7 = 


120.7198 


4869 


B = 


1.2931 


9156 


D = 


-2.0898 


5994 


Br = 


0.8959 


3072 


0,= 


8.0374 


2211 


C = 


0.1377 


8548 


D 2 = 


-7.2966 


5749 


c,= 


14.4323 


0065 


D,= 


1.9772 


8697 


c 2 = 


-79.1011 


8291 


.D 4 = 


-6.3883 


7912 


c,= 


242.3988 


7393 


05= 


6.2341 


2613 



d, mol/1 


P/Pc 


5.0 


0.7373 


6.0 


.8165 


7.0 


.8779 


8.0 


.9237 


9.0 


.9564 


10.0 


.9785 


11.0 


.9918 


12.0 


.9991 


13.0 


1 .0020 


14.0 


1.0021 


14.5 


1.0015 


15.0 


1.0008 


15.5 


1.0001 


16.0 


0.9994 


16.5 


.9989 


17.0 


.9987 


18.0 


.9995 


19.0 


1.0040 


20.0 


1.0149 


21.0 


1.0367 


22.0 


1.0748 


23.0 


1.1367 


24.0 


1.2309 


25.0 


1.3676 


26.0 


1.5584 



Table 4. Comparison ofC* from ft] with calculated values from (5) 






C* at 13.28 mol/1 


T, K 


c° v 








Publ 


Calc 


Diff 


34.216 


12.47 


5.20 


5.30 


-0.10 


36.402 


12.48 


2.79 


3.09 


-.30 


38.384 


12.48 


2.08 


2.21 


-.13 


40.214 


12.49 


1.66 


1.7. 


-.07 


42.107 


12.51 


1.41 


1.41 


-.00 


43.985 


12.52 


1.25 


1.20 


.05 


45.940 


12.55 


1.13 


1.06 


.07 


47.994 


12.59 


1.02 


0.95 


.07 


50.057 


12.63 


0.97 


.88 


.09 


59.868 


13.08 


.80 


.75 


.05 


69.974 


13.98 


.70 


.72 


-.02 


80.206 


15.34 


.69 


.70 


-.01 


90.055 


16.94 


.63 


.68 


-.05 



Note: C* in J/mol K. 
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Table 5. Comparison of C* from [6] with calculated values from (5) 



e «o-(i-i/T) = l _ ( cr _i)8/3. e -o-/(o-)^ 







c\ 


* at 16.18 mo 


/l 


r, k 


c° v 












Publ 


Calc 


Diff 


33.889 


12.47 


6.17 


6.07 


0.10 


36.035 


12.48 


2.73 


3.19 


-.46 


38.194 


12.48 


1.90 


2.13 


-.23 


42.244 


12.51 


1.30 


1.30 


.00 


44.245 


12.53 


1.16 


1.11 


.05 


46.166 


12.55 


1.08 


0.99 


.09 


49.986 


12.63 


0.98 


.87 


.11 


59.894 


13.08 


.84 


.81 


.03 


69.931 


13.98 


.85 


.83 


.02 


80.052 


15.31 


.79 


.83 


-.04 


89.943 


16.94 


.76 


.81 


-.05 



Note: C* in J/mol K. 



Table 6. Comparison ofC* from [2] with calculated values from (5) 







c t 


* at 15.50 mo] 


/l 


T, K 


c; 












Publ 


Calc 


Diff 


33.000 


12.47 


7.24 


13.09 


-5.85 


34.000 


12.47 


5.44 


5.91 


-0.47 


35.000 


12.48 


3.69 


4.23 


-.54 


36.000 


12.48 


2.89 


3.30 


-.41 


37.000 


12.48 


2.16 


2.69 


-.53 


38.000 


12.48 


2.08 


2.26 


-.18 


39.000 


12.49 


1.84 


1.94 


-.09 


40.000 


12.49 


1.65 


1.70 


-.05 


42.000 


12.51 


1.37 


1.36 


.01 


44.000 


12.52 


1.20 


1.15 


.04 


46.000 


12.55 


1.08 


1.02 


.06 


48.000 


12.59 


1.00 


0.93 


.07 


50.000 


12.63 


0.96 


.88 


.08 


55.000 


12.81 


.91 


.81 


.10 


60.000 


13.08 


.84 


.80 


.04 


65.000 


13.47 


.80 


.80 


-.00 


70.000 


13.98 


.77 


.80 


-.03 


75.000 


14.59 


.77 


.80 


-.04 


80.000 


15.31 


.76 


.80 


-.04 


85.000 


16.10 


.77 


.79 


-.02 


90.000 


16.94 


.79 


.78 


.01 


95.000 


17.82 


.81 


.77 


.05 


100.000 


18.69 


.84 


.75 


.09 



Note: C* in J/mol K. 

The term l/x 2 in (7) is qualitatively similar to ^ (p, T) 
at T > 6. Equation (7), however, suffers from coefficient 
D(p) which is negative at all densities [18]. 

The virial equation [18] probably can be related to 
(5) only at such high temperatures that the last two 
terms of (5) are negligible. Figure 2 suggests that 0(p) 
may not be analytic at p = 0. 

7. Appendix. Temperature-density Formula for 
Coexisting Liquid and Vapor 

This relation from [13] is in reduced variables, 
t = T/T c , and a = djd c , 



where 



/(cr) s ai + a2 . (r -i/8 + a8 ., cr i/3' + 

a 4 -cr 8/ ' 3 -f a 5 * o- 3 + a 6 ' cr 7 . 

Values for the constants are, 



a = 2. 91 Ml 


a4 = 


0.620 


5095 


a x = 2.030 0583 


05=- 


-0.421 


81139 


a 2 = - 0.058 7951 


a 6 = 


0.682 


5745 • 10 


a 3 = - 1.856 5706 









The temperature-density relation for coexisting 
liquid and vapor is central to the present equation of 
state. Simplification and increase of accuracy in the 
description of this relation therefore is valuable, as 
obtained with the following form using symbols defined 
above, 

(1/t — l) = |cr — 1| 3 - F(p), 

F(p) =a x - loge (l/p) + a 2 + a 3 P + - • . + a 9 * P 7 - 

Constants used and found for this expression are: 

T c = 32.953 K,d c = 15.59 mol/1, 



d\ = 


0.300 6242 


a 6 = - 


-237.846 


2970 


a 2 = ~ 


- 0.567 0332 


a 7 = 


331.544 


2639 


a 3 = 


3.838 7436 


as = - 


-231.236 


7070 


d4 = - 


-19.378 2878 


a 9 = 


64.190 


7466 


a 5 = 


89.909 2005 









Maximum relative deviations for r are 0.03 percent as 
compared with 0.52 percent in [13]. The overall rms 
relative deviation is 0.016 percent as compared with 
a mean relative deviation of 0.062 percent in [13]. 
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